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Abstract 

The  Goals  algorithm  for  adaptive  modeling  is  extended  to  the  case  of  discrete  mod¬ 
els  such  as  those  characterized  by  lattice  structures  (or  the  static  behavior  molecular 
systems)  by  using  surrogate  models  obtained  from  continuum  approximations  of  the 
lattice.  The  result  is  a  technique  that  could  provide  for  scale-bridging  between  con¬ 
tinuum  models  and  atomistic  models,  although  the  present  development  concerns  only 
simple  algebraic  systems.  An  example  is  provided  in  which  quantities  of  interest  in 
a  large  number  of  degrees  of  freedom  system  are  computed  to  a  preset  tolerance  in 
relatively  few  low-order  approximations. 

Keywords:  multiscale  models,  adaptive  modeling,  error  estimates,  continuum/lattice 
models. 


1  Introduction 


In  earlier  work  [6,  8,  4],  we  developed  techniques  for  assessing  model  and  approximation  error 
in  solid  and  continuum  mechanics,  and  methods  for  adapting  the  mathematical  models  of 
certain  events  to  control  modeling  error.  These  techniques  fall  under  the  general  heading  of 
hierarchical  modeling  and  Goal-oriented  a  posteriori  error  estimation  and  model  adaptivity. 
A  general  theory  for  estimating  modeling  error  was  given  in  [4].  In  the  present  paper, 
we  extend  this  approach  to  the  equilibrium  analysis  of  “atomic”  lattices,  where  the  base 
model  is  defined  by  a  regular  periodic  lattice  and  the  surrogate  models  are  obtained  from 
continuum  models  characterized  by  PDE’s;  precisely  the  opposite  situation  is  encountered 
in  modeling  micro-scale  effects  in  multi-phase  heterogeneous  materials  [6,  8] .  An  important 
byproduct  of  our  analysis  is  a  new  approach  to  modeling  the  transition  from  continuum 
models  to  “molecular”  models.  In  particular,  we  present  techniques  in  which  the  notion  of 
convergence  of  finite  element  models  of  the  continuum  is  unambiguously  preserved,  and  in 
which  convergence  of  the  base  molecular  or  lattice  model  is  trivial.  Thus,  the  notion  of  a 
continuum  atomistic  interface  and  coupling  is  precisely  determined  in  a  rigorous  way. 

The  paper  is  organized  as  follows:  In  Section  2,  we  provide  a  brief  summary  of  the  procedure 
proposed  in  [4]  to  estimate  the  modeling  errors  with  respect  to  some  quantity  of  interest.  We 
present  in  Section  3  a  model  problem  for  static  equilibrium  of  a  material  lattice  and  introduce 
in  Section  4  a  simple  problem  which  will  be  later  referred  to  as  the  continuum  model.  The 
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continuum  model  is  given  in  terms  of  a  partial  differential  equation  and  will  be  considered 
here  as  a  surrogate  problem  of  the  lattice  model  problem.  We  apply  in  Section  5  the 
technique  described  in  Section  2  to  the  continuum/lattice  problems  for  estimating  modeling 
errors,  i.e.  the  difference  between  the  solutions  of  the  continuum  and  lattice  models,  and 
describe  an  adaptive  algorithm  for  controlling  modeling  errors.  We  demonstrate  in  Section  6 
the  efficiency  of  the  adaptive  algorithm  on  two  numerical  examples.  Conclusions  are  given 
in  Section  7. 


2  Modeling  Errors  in  Quantities  of  Interest 


The  abstract  setting  for  the  theory  of  Goal-oriented  estimation  of  modeling  error  is  this:  we 
wish  to  find  a  vector  u  in  a  topological  vector  space  V  such  that 

B(u,v)  =  F(v)  Vv  €  V  (1) 

where  £?(•,•)  is  a  bilinear  (or  semilinear)  form  on  V  and  F  is  a  linear  functional  on  V .  Of 
interest  is  a  quantity  of  interest  Q(u),  Q  being  a  functional  on  V.  To  characterize  Q,  we 
seek  an  influence  vector  p  which  is  a  solution  of  the  dual  problem, 

B(v,p)  =  Q(v)  Vv  €  V  (2) 

Problem  (1)  is  called  the  primal  base  problem;  problem  (2)  the  dual  base  problem.  If  £?(•,  •) 
is  nonlinear  in  u,  we  replace  (2)  by  _B'(u;v,p)  =  Q'(u;v),  where  £?'(■;  •,•)  and  Q'(-\  •)  are 
functional  (Gateaux)  derivatives  of  £?(•,•)  and  Q(-)  respectively. 

The  basic  idea  is  that  (1)  and  (2)  are  intractable:  too  complex  or  too  large  to  be  solved  by 
any  practical  means.  Therefore,  we  seek  simpler  surrogate  problems  that  are  tractable  and 
of  the  form 


B0(u0,v)  =  F(v)  Vv  e  V 
#o(v,po)  =  Q(v)  Vv  G  V 

being  a  new  bilinear  form.  In  [4]  (see  also  [2],  it  is  shown  that 

Q(u)  -  Q(u0)  =  R(u0,  p)  (4) 

where  i?(uo,v)  is  the  residual  functional 

R( u0,  v)  =  -F(v)  -  B( u0,  v)  (5) 

If  the  problem  is  nonlinear,  the  right-hand-side  of  (4)  is  replaced  by  R( uo,  p)  +  r,  where  r  is  a 
functional  involving  quadratic  and  higher  terms  in  the  errors  eo  =  u— Uo  and  eo  =  p  — Po-  In 
applications,  we  often  neglect  such  higher  order  terms,  and  use  i?(uo,  p)  as  an  approximation 
of  Q(u)  —  Q ( Uo ) .  The  use  of  such  relations  to  develop  adaptive  modeling  schemes  for 
estimating  and  controlling  errors  in  modeling  heterogeneous  materials  is  described  in  [6,  8]. 
Our  aim  is  to  use  (4)  as  a  basis  for  estimating  and  controlling  errors  in  models  involving 
a  combination  of  discrete  and  continuum  attributes.  Here,  the  base  model  is  one  provided 
by  the  discrete  lattice  structure,  representing,  for  instance,  a  molecular  dynamics  system  in 
equilibrium,  and  in  which  the  surrogate  model  is  to  be  constructed  from  a  continuum  model 
of  the  material. 


2 


Figure  1:  A  regular  periodic  lattice  in  R2. 


3  A  Model  Lattice  Problem 


To  fix  ideas,  we  consider  a  simple  case  of  static  equilibrium  of  a  material  lattice  with  regular 
periodic  micro-structure.  The  regular  lattice  L  has  lattice  uniform  width  H  and  extends 
indefinitely  in  Rd,  d  =  2  or  3.  For  the  sake  of  simplicity,  we  only  consider  two-dimensional 
domains  here,  i.e.  d  =  2,  as  indicated  in  Fig.  1.  All  results  presented  in  this  report  are 
easily  extended  to  three-dimensional  problems.  The  goal  is  to  determine  a  discrete  scalar 
field  u  which  takes  on  values  u-h3  at  each  lattice  point  x,  j ,  (i,j)  £  Z2,  and  which  satisfies 
the  equilibrium  equation 

X  ahJUi,j  ~  hi  =  0  (6) 

(»,i)6Z2 

where  fij  =  /(xjj),  /  being  an  (i,  j)-periodic  prescribed  force  field,  and  ciij  are  appropriate 
constants. 

To  further  restrict  the  problem  while  also  establishing  conditions  sufficient  to  guarantee  the 
existence  of  solutions  to  (1),  we  shall  make  the  following  simplifications: 

1.  A  subdomain  Q  =  (0,  a)  x  (0,  b )  C  R2  is  identified  with  origin  (0, 0)  at  a  lattice  point, 
as  shown  in  Fig.  2,  and  the  lengths  a  and  b  are  naturally  chosen  as  multiples  of  H. 
Note  that  in  more  general  settings,  O  can  be  a  more  complex  domain  than  just  a 
rectangle  and  may  depend  on  the  structure  of  the  lattice  as  well. 

2.  The  boundary  of  is  denoted  dQ  and  Uij  =  0  at  points  Xjj  £  dfi. 

3.  fij  is  given  at  all  points  x,j  £  O. 

4.  The  coefficients  representing  the  interactions  of  adjacent  ’’atoms”  at  sites  in  the  lattice 
are  given. 

For  simplicity,  one  example  is  the  difference  stencil, 

~Jj2Uid-^  +  ~Jj2Uid  ~  —  ~  fi,j  =  (7) 

which  we  write  in  matrix  form  as 

Au  =  f  (8) 
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Our  goal  is  not  to  calculate  u,  but  some  functional  of  u  defined  by  a  quantity  of  interest 
Q{  u).  For  example,  if  w  is  a  patch  of  four  cells  surrounding  a  particular  lattice  point  x* 
(see  Fig.  2),  we  may  be  interested  in  the  weighted  average  value  of  u  over  the  cell: 

Q( u)  =  (wfei  +  Uk2  +  Uk3  +  Uk4)  +  2(uk5  +  Uk6  +  Uh7  +  Uk3)  +  4ufeg 

9  (12) 

=  'y  ]  wkiUki 

i=i 

Wk,  being  the  indicated  weights.  This  corresponds  to  the  average  of  the  interpolant  of  u 
over  the  patch  u>.  The  problem  of  finding  a  vector  p  S  14  such  that 

£(v,  p)  =  Q(v)  Vv  G  VL  (13) 

is  then  the  base  dual  problem  for  the  functional  Q,  and  p  is  the  dual  solution  or  influence 
vector  for  the  lattice  problem. 
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4  A  Continuum  Surrogate  Problem 


We  propose  here  a  simple  continuum  model  deduced  from  the  lattice  model.  Our  principal 
motivation  in  this  paper  is  to  focus  on  the  adaptive  algorithm  rather  than  modeling  issues. 
Various  techniques  have  been  suggested  to  derive  continuum  models,  see  for  example  [3,  7]. 
For  the  lattice  model  (8),  in  the  limit  as  H  — >  0,  we  consider  here  the  following  continuum 
model  of  the  problem: 


Find  uq  (E  V  such  that  B0(uq,v)  =  F0(v)  Vu  €  V 

(14) 

and 

Find  po  €  V  such  that  B0(v,po)  =  Qo(v)  Vu  €  V 

(15) 

where  now 

Bq(u,v)=  /  S7u-S7vdxdy 

Jn 

(16) 

F0(v)  =  [  fv  dxdy 

Jci 

(17) 

Qo(v)  =  i— r  /  v  dxdy 

M  L 

(18) 

V  =  {v  G  B1(fl)  :  v  =  0  on  =  Hq(Q) 

(19) 

Clearly,  in  the  limit  as  H  — >  0,  the  lattice  model  (8)  converges  to  the  model  of  the  Poisson 
problem  —Art  =  /  in  O,  u  =  0  on  dil.  Of  course,  the  spaces  Vl  and  V  are  now  incompatible, 
Vl  being  of  dimension  N  and  V  infinite  dimensional. 

We  must  next  map  the  continuum  description  onto  the  lattice  model  to  be  able  to  compare 
errors  in  the  quantities  of  interest.  This  is  done  as  follows.  Let  II  :  V  — >  Vl  denote  a 
(collocation)  mapping  from  V  onto  Vl  defined  by 

nw  =  v  =  {vfc}f=1 

Vk=Wp(xk)=  /  fep(x,xfc)w(x)  dx 

J  Bp  (xfc)  (20) 

p<^H 

Xfc  =  lattice  site  in  L 

Here  _Bp(xfc)  is  a  ball  of  radius  p  centered  at  the  lattice  node  Xfc,  fcp(x,  Xfc)  is  a  smooth  kernel 
vanishing  outside  _Bp(xfc)  and  normalized  so  that 

/  fcp(x,xfc)  dx  =  1. 

J  Bp(xfc) 

Thus,  wp  is  the  mollifier  (or  mollification)  of  w.  The  operator  n  thus  produces  an  N -vector 
in  Vl  defined  on  the  lattice  whose  components  are  simply  the  values  of  the  functions  w  in 
V,  smoothed  so  that  their  pointwise  values  are  well-defined  (the  functions  w  €  V,  being  in 
H 1(f2),  are  not  necessarily  continuous  for  2D  and  3D  domains).  In  practice,  we  replace  wp 
by  finite  element  approximations  w h  for  mesh  size  h  sufficiently  small.  This  is  discussed  in 
Section  5,  to  follow. 

Applying  the  ideas  of  Section  2,  we  now  derive  an  estimate  for  the  modeling  error.  Let 
mq  G  V  be  the  unique  solution  of  (14)  and  p  €  14  be  the  influence  vector  for  Q  satisfying 
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(13).  Then,  according  to  (4),  the  modeling  error  eo  =  u  —  ILuo  =  u  —  Uo  in  the  quantity  of 
interest  Q  is 

<2(eo)  =  Q(u)  -  Q(u0)  =  i?(u0,  p).  (21) 

We  can  rewrite  the  residual  as: 

#(u0,p)  =  i?(u0,p0)  +  f?(u0,£o)  (22) 

=  i?(u0,p0)  +  B(e0,EQ)  (23) 

where  po  =  Ilpg,  Po  being  the  influence  function  for  the  surrogate  problem  (15),  and  eg  = 
p  —  p0.  Note  that  (23)  follows  from  (22)  by  making  use  of  the  error  equation  B(e0,v)  = 
R( Uo,  v),  Vv  G  V  and  by  taking  v  =  e0.  An  estimate  of  the  modeling  error  is  then  obtained 
by  neglecting  the  higher  order  term  B(e0,e0)  so  that  Q(e0)  ns  i?(u0,p0),  or  by  computing 
bounds  on  B(e0,  £o)  as  in  [6],  or  computing  the  exact  dual  solution  vector  p. 


5  Approximations  of  the  Continuum  Model 

5.1  Discrete  approximations 

In  general,  the  continuum  models  (14)  and  (15)  cannot  be  solved  exactly,  and  numerical 
approximations,  such  as  finite  element  approximations,  of  uq  and  po,  must  be  obtained. 
Then,  instead  of  (20),  a  more  direct  construction  of  a  mapping  from  V  onto  Vl  can  be 
defined. 

Let  {Ph}  denote  a  family  of  partitions  of  fi  into  finite  elements  K  each  being  the  image 
under  invertible  maps  Fk  of  a  master  element  I\  over  which  polynomial  test  functions  are 
defined  (see,  e.g.  [1]).  We  shall  assume  that  the  partitions  are  regular  and  that  the  usual 
interpolation  properties  of  finite  element  interpolation  of  functions  in  Sobolev  spaces  are 
in  force  [5].  In  this  way,  we  generate  a  family  of  finite-dimensional  subspaces  {Vh}  of  V, 
with  IJ/i>o  Vh  everywhere  dense  in,  say,  H1( fi).  For  a  given  subspace  Vh,  the  finite  element 
approximations  of  (14)  and  (15)  are  thus 

B0(u%,vh)  =  F0(vh)  \/vh  G  Vh 
B0(vh,p%)  =  Q0(vh)  Vvh  G  y'1 

Instead  of  (20),  we  introduce  the  “collocation”  operator  n^:  Vh  —>  Vl  such  that 

Uhvh  =  vh  =  {(vh)k}^i 

(vh)k  =  Vh(xk)  (26) 

Xfc  =  lattice  site  in  O 


(24) 

(25) 


5.2  Estimate  of  modeling  error 

To  estimate  the  error  in  Q(uq),  we  proceed  as  before.  Using  Uq  =  IEwq  and  pj  =  n^pg 
instead  of  u0  and  po,  respectively,  the  estimate  (21)  becomes 

Q(u)  —  Q(ug)  =  i?(u^,p) 

=  -^(U0)Po)  +  —  Ug,p  —  Pg)  (27) 
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(28) 


Neglecting  the  higher  order  term  in  u  ufj,  we  then  obtain  the  error  estimate: 

£  =  #K,pft)«Q(u)-Q(<) 

It  is  important  to  note  that  there  is  no  connection  between  the  lattice  dimension  H  and 
the  mesh  size  h.  Convergence  of  the  finite  element  approximation  is  unambiguous:  as 
h  — >  0,  Uq  — >  uo  in  H1(  12),  independently  of  H  (for  quasi-uniform  refinements,  with  h  = 
maxftx,  hx  =  dia(AT)). 


5.3  Adaptive  modeling  strategy 


We  now  describe  an  adaptive  modeling  strategy  aimed  at  reducing  computational  costs. 
For  instance,  let  us  suppose  that  we  cannot  solve  the  iV-order  system  (11),  N  being  a  large 
number,  but  that  we  wish  to  compute  for  the  vector  u  the  quantity  of  interest  Q(u)  to 
within  a  tolerance  7 toi,  i.e. 

IQ(u)  ~Q(uq)I 

|Q(u)|  “  toh 

We  propose  a  Goal-oriented  hierarchical  modeling  algorithm  that  employs  concepts  similar 
to  those  developed  in  [6,  8]  for  computing  local  features  in  a  lattice.  Again,  let  L  denote 
the  lattice  of  N  sites  covering  the  domain  O  with  lattice  width  H.  In  addition,  let  fi  be 
partitioned  into  K  cells  of  size  larger  than  H  (see  Fig.  3).  The  error  estimate  £  can  then 
be  decomposed  into  contributions  over  each  cell  as 

K 

£  =  ft(uo,Po)  =  X^(uO’Po)  (29) 

n—  1 


where  1Zn  ( Uq  ,  p(j  )  denotes  the  contribution  in  the  nth  cell  to  the  total  error  £.  More  explic¬ 
itly,  recall  that 

^(uo,Po)  =  ^(Po)-^(uo,Po) 

=  PoTf  —  PoTj4uq 
=  p&r(f-Au§) 

N  /  N  \ 

=  po k  ( ffc  -  Akiuoi ) 

k= 1  V  1=1  ) 


Then  the  contribution  Sn  from  the  nttl  cell  0n  is  computed  according  to 

N 

en  =  Kn(u!;,pS)  =  Z2/% 


Po k  ffc  -  >  Afc;u; 


N 

'  '  ‘  ..h 

klu0l 


Z _ j 


(30) 


fc=i 


1=1 


where  the  coefficients  /3£  depend  on  the  location  of  the  sites  with  respect  to  the  cell  0„. 
Let  6k  be  the  number  of  cells  that  contain  point  x*,.  This  number  can  take  on  the  values 
4,  2,  or  1,  depending  on  whether  xj,  is  an  interior  site  of  L,  or  lies  on  one  of  the  boundary 
edges  of  L,  or  is  one  of  the  four  corners  of  L ,  respectively.  Then 


Pk 


0  if  xk  0„ 
9 (T1  if  xk  G  0„ 


The  algorithm  for  adaptively  modeling  the  lattice  model  with  a  continuum  model  is  as 
follows: 
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1.  Set  s  =  0 


2.  Replace  the  lattice  model  with  the  continuum  model  (14)  and  solve  for  a  (very  accu¬ 
rate)  finite  element  approximation  u g  of  Uq.  Compute  also  an  approximation  p ft  of 
Po- 

3.  Evaluate  Uq  =  n^Ug,  pg  =  HhPo,  and  the  quantity  of  interest  Q(ug). 

4.  Estimate  the  error  in  Q(ug  )  using  (28),  i.e.  compute  =  i?(ug,pg). 

5.  Check  whether  |£(°) |/|Q(ug)  +  £®|  <  7 toil  If  yes,  stop.  If  not,  start  the  adaptive 
process. 

6.  Set  s  =  s  +  1 

(a)  Compute  the  cell  contributions  £n  =  Un(uJ_i, Ps_i),  n  =  1  as  in  (30). 

Construct  the  patch  Ls  made  of  the  cells  that  have  not  been  already  refined  and 
that  satisfy 

16,1  . 

- FFT  -  a 

max„  \  tn\ 

where  a  is  a  user-defined  parameter. 

(b)  Solve  the  reduced  lattice  problem  on  Ls  for  lattice  vectors  uj  and  p(b  with 
uj  =  Ug  and  p((  =  pg  at  the  points  outside  of  Ls. 

(c)  Estimate  the  error  in  Q(uJ),  i.e.  compute  £^s>  =  R( u(f,Pg)  . 

(d)  Check  whether  |£(s)  |/|Q(Ug)+£W  |  <  7  to;?  If  the  answer  is  negative,  go  to  step  6. 
If  the  answer  is  affirmative,  stop  the  adaptive  process. 

In  this  manner,  we  create  a  sequence  of  surrogate  problems  Vs,  whose  solutions  are  given  by 
(uj ,  Pg),  s  =  0,  1,  ....  The  initial  surrogate  problem  Vo  is  obtained  using  the  continuum 
model  only.  The  subsequent  problems  Vs,  s  =  1,  2,  . . . ,  combine  the  solutions  from  the 
continuum  model  and  from  the  lattice  model  used  only  in  the  reduced  lattice  Ls.  Obviously, 
this  is  only  one  example  of  many  possible  variants  of  the  adaptive  process.  Note  that  the 
overall  approach  is  also  reminiscent  of  the  multigrid  method  in  which  the  lattice  model 
stands  for  the  the  fine-scale  grid,  and  the  surrogate  model  for  the  coarse-scale  grid  model. 


6  Numerical  Examples 


To  demonstrate  the  Goals  algorithm  for  adaptive  modeling,  we  consider  examples  of  regular 
lattices  L  on  the  square  domain  £1  =  (0,  l)2  with  lattice  width  H  =  (m  —  1)_1,  m  =  11,  21, 
. . . ,  61.  The  lattice  L  is  then  made  of  N  =  m2  sites  such  that: 

L  =  {(27,2/7):  bJ  e  N,  0  <  i,  j  <  TO  -  1,  Xi  =  iH:  yj  =  jH} 

and  will  be  referred  to  as  a  m  x  m  lattice.  We  emphasize  that  the  lattice  solution  is  our 
solution  of  reference,  although  in  these  examples,  the  scales  provided  by  the  finite  element 
solutions  are  finer.  Again  our  only  concern  here  is  to  demonstrate  the  feasibility  of  the 
adaptive  strategy. 

In  all  examples,  we  partition  the  domain  f2  into  25  cells  of  dimension  0.2  x  0.2  (see  Fig.  3) 
and  the  quantity  of  interest  is  the  weighted  average  over  the  square  domain  u>  of  dimension 
2  H  x  2 H  located  at  the  center  of  the  lattice,  as  defined  in  (12).  In  other  words,  the 
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Figure  3:  A  lattice  in  R2. 

quantity  of  interest  is  the  weighted  average  of  u  at  the  lattice  points  Xjj  where  i,j  £ 
{(to  —  l)/2  —  1,  (to  —  l)/2,  (to  —  l)/2  +  1}-  Note  that  the  quantity  of  interest  depends 
indirectly  on  the  lattice  width  H . 

The  exact  solution  to  the  continuum  model  problem  is  given  by 

—Auo  =  /,  in  no  =  0  on  dil 

where  the  datum  /  is  determined  such  that  Uo(x,  y)  is  known  analytically.  We  will  consider  in 
the  following  two  different  functions  uq.  The  solution  of  the  Poisson  problem  is  approximated 
using  a  mesh  of  200  x  200  bilinear  elements.  In  all  numerical  experiments,  the  parameter  a 
introduced  in  step  6(a)  of  the  adaptive  algorithm  will  be  chosen  to  be  0.5. 


6.1  Example  one 

In  this  example,  the  continuum  and  lattice  problems  are  set  up  in  such  a  way  that  the  exact 
continuum  solution  is  given  by 

u0(x,y)  =  256x2(1  -  x)2y2(l  -  yfe"100^-0-25)^-0-25)2) 

The  function  u o  is  essentially  a  differenced  Gaussian  centered  at  the  point  (0.25,  0.25)  as 
shown  in  Fig.  4.  The  term  a;2(l  —  x)2y2(  1  —  y)2  ensures  that  uq  and  the  derivatives  of  uq 
are  zero  on  the  boundary  of  O.  The  corresponding  lattice  solutions  are  shown  in  Fig.  5 
for  the  11  x  11  and  31  x  31  lattices.  We  show  in  Figs.  6  and  7  the  continuum  influence 
function  po  and  lattice  influence  influence  p  associated  with  the  weighted  average  defined 
above,  for  the  11  x  11  and  31  x  31  lattices,  respectively.  As  expected,  we  observe  that  the 
influence  function  becomes  more  localized  as  the  lattice  width  decreases  since,  in  all  cases, 
the  quantity  of  interest  is  defined  with  respect  to  the  nine  lattice  points  at  the  center  of  the 
lattice. 
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Figure  4:  Accurate  approximation  of  the  continuum  solution  u o  for  example  1. 


Figure  5:  Lattice  solutions  of  example  1  computed  on  the  11  x  11  and  31  x  31  lattices. 


In  the  next  set  of  numerical  experiments,  we  test  the  adaptive  modeling  strategy  proposed 
in  Section  5.3.  We  set  the  tolerance  on  the  error  7 toi  to  5  percent.  Figures  8  and  9  show 
the  sequence  of  refinements  for  the  21  x  21  lattice.  In  these  plots,  and  in  all  similar  plots 
that  follow,  the  grey  area  represents  the  subregion  in  which  the  solution  has  been  computed 
using  the  continuum  model,  that  is,  wo-  The  complementary  subregion  displaying  the  lattice 
corresponds  to  the  reduced  lattice  Ls ,  where  s  is  the  number  of  iterations  already  performed, 
as  explained  in  Section  5.3.  The  number  l  <  s  in  the  various  cells  indicates  the  iteration  at 
which  the  cell  has  been  added  to  the  lattice  Ls.  In  this  case,  eight  refinements  were  necessary 
to  achieve  a  tolerance  of  5  percent  which  resulted  in  all  but  three  cells  being  refined.  Note 
that  all  plots  are  symmetric  with  respect  to  the  diagonal  x  =  y,  as  expected. 

We  show  in  Fig.  10  the  final  configurations  obtained  by  adaptive  modeling  for  all  lattices. 
We  observe  that  eight  iterations  are  needed  for  the  11  x  11  and  21  x  21  lattices  to  achieve  a 
tolerance  of  5  percent  on  the  error  in  the  quantity  of  interest.  For  all  the  other  lattices,  the 
algorithm  yields  the  same  final  configuration  after  5  iterations,  but  note  that  the  cells  are 
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Figure  6:  Influence  functions  associated  with  the  11  x  11  lattice:  (left)  continuum  solution 
Po  and  (right)  lattice  solution  p. 


Figure  7:  Influence  functions  associated  with  the  31  x  31  lattice:  (left)  continuum  solution 
po  and  (right)  lattice  solution  p. 


not  necessarily  refined  in  the  same  order  in  all  cases. 

In  Fig.  11,  we  show  the  effectivity  indices  and  the  evolution  of  the  relative  errors  with  respect 
to  the  number  of  refinements  obtained  by  adaptive  modeling.  The  effectivity  indices  provide 
a  measure  of  the  quality  of  the  error  estimates  and  are  calculated  as  the  ratio  of  the  error 
estimate  £ ^  at  iteration  s  and  the  exact  error,  i.e. 

|g(s)l 

|Q(u)-Q(uJ)| 

Thus,  an  effectivity  index  of  exactly  one  indicates  that  the  error  estimator  perfectly  estimates 
the  exact  error.  From  Fig.  11,  we  observe  that  the  effectivity  indices  oscillates  between  0.95 
and  1.1  except  for  one  value.  The  rates  of  convergence  are  mostly  monotonic  apart  from  the 
case  of  the  11  x  11  lattice.  In  this  case,  we  believe  that  some  internal  boundary  conditions 
obtained  from  the  continuum  solution  are  in  error  and  greatly  increase  the  relative  error  in 
the  quantity  of  interest  until  these  boundaries  become  part  of  the  lattice  by  including  the 
adjacent  cells. 

Finally,  we  show  in  Fig.  12,  for  illustration  purposes,  the  solutions  obtained  on  the  reduced 
lattices  L 3  (after  the  third  iteration)  for  the  11  x  11  and  31  x  31  lattices.  These  plots  can 
be  qualitatively  compared  to  the  solutions  that  were  computed  on  either  the  full  continuum 
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Figure  8:  Sequence  of  surrogate  models  obtained  by  adaptive  refinement  for  the  21  x  21 
lattice  (iterations  s  =  1  to  4).  In  these  plots,  the  grey  area  represents  the  subregion  in 
which  the  solution  is  obtained  using  the  continuum  model.  The  complementary  subregion 
corresponds  to  the  reduced  lattice  Ls.  The  number  l  <  s  in  the  cells  indicates  the  iteration 
at  which  the  cell  has  been  added  to  the  lattice  Ls. 


model  or  the  full  lattice  model;  see  Figs.  4  and  5. 


6.2  Example  two 

We  repeat  in  this  section  the  same  types  of  numerical  experiments  as  before  but  for  a  different 
exact  solution  of  the  continuum  model.  Here  the  solution  uq  (see  Fig.  13)  is  obtained  by 
superposition  of  two  “Gaussian  functions”  centered  at  the  points  (0.25, 0.25)  and  (0.50,  0.80) 
such  that: 


u0(x,y)  =128cc2(1  -  a)V(l  -  3,)2e-ioo((*-o.25)3+(,,-o.25)a) 

+  512*2(1  -  X)2V2{  1  -  2/)2e-100((x-0.50)W(y-0.80)2) 

In  this  case,  the  weighted  average  error,  still  computed  over  the  nine  points  at  the  center  of 
the  lattices,  is  essentially  influenced  by  the  sources  of  error  generated  at  the  location  of  the 
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Figure  9:  Sequence  of  surrogate  models  obtained  by  adaptive  refinement  for  the  21  x  21 
lattice  (iterations  s  =  5  to  8).  In  these  plots,  the  grey  area  represents  the  subregion  in 
which  the  solution  is  obtained  using  the  continuum  model.  The  complementary  subregion 
corresponds  to  the  reduced  lattice  Ls.  The  number  l  <  s  in  the  cells  indicates  the  iteration 
at  which  the  cell  has  been  added  to  the  lattice  Ls. 


two  Gaussian  functions.  The  corresponding  lattice  solutions  are  shown  in  Fig.  14  for  the 
11  X  11  and  31  x  31  lattices.  Since  the  quantities  of  interest  are  the  same  as  in  example  1, 
so  are  the  associated  influence  functions. 

For  testing  the  adaptive  modeling  strategy  proposed  in  Section  5.3,  we  select  two  tolerances, 
i.e.  7 toi  =  10  percent  and  7 toi  =  5  percent.  Fig.  15  shows  the  sequence  of  refinements  for  the 
21  x  21  lattice.  The  algorithm  automatically  performed  five  refinements  in  order  to  achieve 
a  tolerance  of  10  percent.  We  observe  that  the  procedure  essentially  refines  the  cells  that 
cover  the  peaks  of  the  Gaussian  functions  and  the  features  of  the  quantity  of  interest.  The 
final  configurations  for  all  lattices  are  shown  in  Fig.  16  for  7 toi  =  10  percent  and  in  Fig.  17 
for  7 toi  =  5  percent.  We  observe  that,  apart  from  the  11  x  11  lattice,  two  to  three  extra 
iterations  are  necessary  to  reduce  the  error  from  10  percent  to  5  percent.  The  evolution 
of  the  effectivity  indices  and  relative  errors  with  respect  to  the  number  of  refinements  are 
displayed  in  Figs.  18  and  19.  The  convergence  results  suggest  that  the  adaptive  algorithm 
could  be  improved  for  our  particular  choice  of  quantity  of  interest.  The  erratic  behavior  prior 
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to  a  dramatic  reduction  in  error  may  be  attributed  to  the  fact  that  the  sequence  of  model 
enrichments  employs  local  error  measures  rather  than  global  error  measures.  Regardless  of 
the  rate  of  convergence,  quite  acceptable,  indeed  excellent,  estimates  of  the  modeling  error 
are  obtained.  A  more  detailed  study  of  various  adaptive  strategies  and  convergence  will  be 
undertaken  for  more  complex  systems  in  future  studies.  Finally  the  solutions  obtained  on 
the  reduced  lattices  L2  for  the  11  x  11  and  31  x  31  lattices  are  given  in  Fig.  20  and  should  be 
compared  with  the  solutions  of  Figs.  13  and  14.  These  pictures  clearly  demonstrate  that  the 
proposed  adaptive  procedure  automatically  refines  the  regions  in  which  the  large  sources  of 
modeling  error  contribute  the  most  to  the  error  in  the  quantity  of  interest. 


7  Concluding  Remarks 


We  have  described  in  this  paper  a  method  to  extend  the  Goals  algorithm  for  adaptive  mod¬ 
eling  to  the  case  of  discrete  models  by  using  a  continuum  model  for  the  surrogate  problem. 
The  method  involves  the  derivation  of  a  Goal-oriented  error  estimator  to  obtain  computable 
error  measures  of  local  quantities  of  interest.  The  proposed  adaptive  algorithm  is  an  it¬ 
erative  procedure  which  allows  one  to  determine  the  regions  of  the  domain  in  which  it  is 
necessary  to  use  the  lattice  base  model  to  control  the  error  in  the  quantity  of  interest  to 
within  some  preset  tolerance.  The  solution  of  the  continuum  model  is  merely  used  to  pre¬ 
scribe  the  internal  boundary  conditions  for  the  reduced  lattice  problems.  The  methodology 
was  tested  here  on  two  numerical  examples  and  the  results  clearly  demonstrate  the  great 
potential  of  this  approach.  In  these  examples,  the  continuum  model  is  defined  as  a  Poisson 
problem  and  the  lattice  model  was  derived  by  using  a  five-point  central  difference  stencil 
for  the  Poisson  equation.  The  main  conclusion  from  this  study  is  that  it  is  actually  feasible 
to  automatically  select  the  sites  of  the  lattice  that  need  be  included  in  the  reduced  lattice 
problem  to  obtain  accurate  quantities  of  interest.  Whenever  the  lattice  problem  is  more 
expensive  to  solve  than  the  continuum  problem,  the  use  of  this  approach  would  allow  for 
substantial  cost  reductions.  Our  objective  in  the  near  future  is  to  extend  this  methodol¬ 
ogy  to  molecular  statics  or  molecular  dynamics  for  problems  in  nanomechanics.  The  main 
issue  will  be  to  construct  adequate  surrogate  problems  from  the  lattice  problems  for  multi¬ 
scale  modeling  and  to  deal  with  the  loss  of  scale  information  due  to  the  transition  from  the 
molecular  model  to  the  continuum  model. 
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Figure  10:  Adaptive  modeling  for  the  latt: 
For  each  case  the  final  configuration  whicl 
the  quantity  of  interest  is  shown. 


Figure  11:  Effectivity  indices  (left)  and  relative  errors  (right)  versus  the  number  of  refine¬ 
ments  for  example  1. 


0.4  . 
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Figure  12:  Solutions  obtained  on  the  reduced  lattices  Lj,  (after  iteration  3)  for  the  11  x  11 
lattice  (top)  and  the  31  x  31  lattice  (bottom). 
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Figure  13:  Accurate  approximation  of  the  continuum  solution  ug  for  example  2. 


Figure  14:  Lattice  solutions  of  example  2  computed  on  the  11  x  11  and  31  x  31  lattices. 
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Figure  15:  Sequence  of  surrogate  models  obtained  by  adaptive  refinement  for  the  21  x  21 
lattice  (iterations  s  =  1  to  5).  Again,  the  grey  area  represents  the  subregion  in  which  the 
solution  is  obtained  using  the  continuum  model.  The  complementary  subregion  corresponds 
to  the  reduced  lattice  Ls.  The  number  l  <  s  in  the  cells  indicates  the  iteration  at  which  the 
cell  has  been  added  to  the  lattice  Ls. 
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Figure  16:  Adaptive  modeling  for  the  lattices  11  x  11,  21  x  21,  . . . ,  61x61  for  example  2. 
For  each  case  the  final  configuration  which  achieves  a  tolerance  of  10  percent  in  the  error  in 
the  quantity  of  interest  is  shown. 
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Figure  17:  Adaptive  modeling  for  the  lattices  11  x  11,  21  x  21,  . . . ,  61  x  61  for  example  2. 
For  each  case  the  final  configuration  which  achieves  a  tolerance  of  5  percent  in  the  error  in 
the  quantity  of  interest  is  shown. 
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Figure  18:  Effectivity  indices  (left)  and  relative  errors  (right)  versus  the  number  of  refine¬ 
ments  for  example  2  and  7toi  =  10  percent. 


Figure  19:  Effectivity  indices  (left)  and  relative  errors  (right)  versus  the  number  of  refine¬ 
ments  for  example  2  and  7toi  =  5  percent. 
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Figure  20:  Solutions  obtained  on  the  reduced  lattices  L2  (after  iteration  2)  for  the  11  x  11 
lattice  (top)  and  the  31  x  31  lattice  (bottom). 
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